#include <pqxx/pqxx>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <dfsd.h>
#include <vector>
#include <mfhdf.h>
#include <hdf.h>
#include <gdal.h>

#include <sys/types.h>
#include "libpq-fe.h"

#include <iostream>
#include <sstream>
#include <fstream>

using namespace std;
using namespace pqxx;


typedef struct
{
    float32 lat;
    float32 lon;

    // surface reclectance for bands
    int16 sr_b01; //  Surface_reflectance_for_band_1
    int16 sr_b02; //  Surface_reflectance_for_band_2
    int16 sr_b03; //  Surface_reflectance_for_band_3
    int16 sr_b04; //  Surface_reflectance_for_band_4
    int16 sr_b05; //  Surface_reflectance_for_band_5
    int16 sr_b06; //  Surface_reflectance_for_band_6
    int16 sr_b07; //  Surface_reflectance_for_band_7

    // quality control
    uint32 sr_qc_500m;        //  Surface_reflectance_500m_quality_control_flags
    int16 sr_szen;            //  Solar_zenith
    int16 sr_vzen;            //  View_zenith
    int16 sr_raz;             //  Relative_azimuth
    uint16 sr_state_500m;     //  Surface_reflectance_500m_state_flags
    uint16 sr_day_of_year;    //  Surface_reflectance_day_of_year


} CombiData;

GByte getHDFMODPixelValueT(char *fMODISName, char *bName, int bNum, VOIDP outRaster);
void    findDimension(char *nameInputFile, int **arrayDimension, char *nameBand);
void    getFactor(char * inputFileName, char * band, char *attribute, float **valueToSearch);


int main()//int argc)//, char* argv[])
{

    int tileID = 51027006;
    // input file
    int left_coo = 2232;
    int right_coo = 2394;
    int top_coo = 2225;
    int bot_coo = 2334;
    
    char fMOD09A1Name[100];
    int *dimDataset, dx, dy, k, j, pos;

    GByte resHdf;
    CombiData data;

    // file name
    ifstream ifs;
    char filename[50];
    ifs.open ("/media/MEDIA/DataNCKH/listfile.txt", std::ifstream::in);

    string query;
    stringstream tmp_query;    

    try
    {
        connection C("dbname=mymodisdb user=demo1usr password=demo1pwd hostaddr=127.0.0.1 port=5432");
        if (C.is_open()) {
            cout << "Opened database successfully: " << C.dbname() << endl;
        } else {
            cout << "Can't open database" << endl;
            return -1;
        }
        /*
        // points - table convert x,y in pixel to lat long in wgs84/sinusoidal projection
        for (int xcoo = left_coo ; xcoo <= right_coo ; xcoo++)
        {
            for (int ycoo = top_coo ; ycoo <= bot_coo ; ycoo++)
            {
                /* Create query statement
                tmp_query.str("");
                tmp_query << "INSERT INTO points (px_x, px_y, lat, lon ) VALUES (" << xcoo << "," << ycoo << ", " << 0 << ", " << 0 << ");";
                query = tmp_query.str();
                work W(C);
                W.exec( query );
                W.commit();
            }
        }
        */
        cout << "Create table point successfully. " << C.dbname() << endl;
        
        // loop images : 4 year and 46 images each
        for (int year = 2012; year <= 2012; year++)
        {
            for (int img = 217; img <= 365; img=img+8)
            {

                int16 *sr_b01 = NULL;
                int16 *sr_b02 = NULL;
                int16 *sr_b03 = NULL;
                int16 *sr_b04 = NULL;
                int16 *sr_b05 = NULL;
                int16 *sr_b06 = NULL;
                int16 *sr_b07 = NULL;

                uint32 *sr_qc_500m;
                int16 *sr_szen;
                int16 *sr_vzen;
                int16 *sr_raz;
                uint16 *sr_state_500m;
                uint16 *sr_day_of_year;

                // file name
                ifs.getline(filename, 50);
                stringstream tmp_filename;
                tmp_filename << "/media/MEDIA/DataNCKH/" << year << "/" << filename; 
                string strtmp = tmp_filename.str();
                strcpy(fMOD09A1Name, strtmp.c_str());
                string date_taken = strtmp.substr(40,3);

                // image dimension
                findDimension(fMOD09A1Name, &dimDataset, "sur_refl_b01");
                dy =    (int) dimDataset[0];
                dx =    (int) dimDataset[1];


                // get scale factor
                //float *scale_factor_b01;
                //float *scale_factor_b02;
                //float *scale_factor_b03;
                //float *scale_factor_b04;
                //float *scale_factor_b05;
                //float *scale_factor_b06;
                //float *scale_factor_b07;
                //getFactor(fMOD09A1Name, "sur_refl_b01", "scale_factor", &scale_factor_b01);
                //getFactor(fMOD09A1Name, "sur_refl_b02", "scale_factor", &scale_factor_b02);
                //getFactor(fMOD09A1Name, "sur_refl_b03", "scale_factor", &scale_factor_b03);
                //getFactor(fMOD09A1Name, "sur_refl_b04", "scale_factor", &scale_factor_b04);
                //getFactor(fMOD09A1Name, "sur_refl_b05", "scale_factor", &scale_factor_b05);
                //getFactor(fMOD09A1Name, "sur_refl_b06", "scale_factor", &scale_factor_b06);
                //getFactor(fMOD09A1Name, "sur_refl_b07", "scale_factor", &scale_factor_b07);

                float scale_factor_b01 = 0.0001;
                float scale_factor_b02 = 0.0001;
                float scale_factor_b03 = 0.0001;
                float scale_factor_b04 = 0.0001;
                float scale_factor_b05 = 0.0001;
                float scale_factor_b06 = 0.0001;
                float scale_factor_b07 = 0.0001;

                // get offset
                float *add_offset_b01;
                float *add_offset_b02;
                float *add_offset_b03;
                float *add_offset_b04;
                float *add_offset_b05;
                float *add_offset_b06;
                float *add_offset_b07;
                getFactor(fMOD09A1Name, "sur_refl_b01", "add_offset", &add_offset_b01);
                getFactor(fMOD09A1Name, "sur_refl_b02", "add_offset", &add_offset_b02);
                getFactor(fMOD09A1Name, "sur_refl_b03", "add_offset", &add_offset_b03);
                getFactor(fMOD09A1Name, "sur_refl_b04", "add_offset", &add_offset_b04);
                getFactor(fMOD09A1Name, "sur_refl_b05", "add_offset", &add_offset_b05);
                getFactor(fMOD09A1Name, "sur_refl_b06", "add_offset", &add_offset_b06);
                getFactor(fMOD09A1Name, "sur_refl_b07", "add_offset", &add_offset_b07);

                tmp_query.str("");
                tmp_query << "INSERT INTO images ( tile_id, date_tkn, b01_scale_factor, b01_add_offset, b02_scale_factor, b02_add_offset, b03_scale_factor, b03_add_offset, b04_scale_factor, b04_add_offset, b05_scale_factor, b05_add_offset, b06_scale_factor, b06_add_offset, b07_scale_factor, b07_add_offset ) VALUES (" << tileID << ",'" << year << "." << date_taken << "'," << scale_factor_b01 << "," << *add_offset_b01 << "," << scale_factor_b02 << "," << *add_offset_b02 << "," << scale_factor_b03 << "," << *add_offset_b03 << "," << scale_factor_b04 << "," << *add_offset_b04 << "," << scale_factor_b05 << "," << *add_offset_b05 << "," << scale_factor_b06 << "," << *add_offset_b06 << "," << scale_factor_b07 << "," << *add_offset_b07 << ");"; 
                query = tmp_query.str();
                cout << filename << " processed!" << endl;

                work W(C);
                W.exec(query);
                W.commit();

                // MOD09A1 data
                sr_b01  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b02  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b03  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b04  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b05  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b06  = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_b07  = (int16 *)malloc(sizeof(int16) * dx * dy);

                if (sr_b01 == NULL || sr_b02 == NULL || sr_b03 == NULL || sr_b04 == NULL ||
                    sr_b05 == NULL || sr_b06 == NULL || sr_b07 == NULL){
                    printf ("Malloc Error");
                    return (1);
                }

                // MOD09A1 quality control
                sr_qc_500m  = (uint32 *)malloc(sizeof(uint32) * dx * dy);
                sr_szen     = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_vzen     = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_raz      = (int16 *)malloc(sizeof(int16) * dx * dy);
                sr_state_500m   = (uint16 *)malloc(sizeof(uint16) * dx * dy);
                sr_day_of_year  = (uint16 *)malloc(sizeof(uint16) * dx * dy);

                if (sr_qc_500m == NULL || sr_szen == NULL || sr_vzen == NULL ||
                    sr_raz == NULL || sr_state_500m == NULL || sr_day_of_year == NULL){
                    printf ("Malloc Error");
                    return (1);
                }

                // Get surface relfectance at 7 bands
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b01", 0, sr_b01);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b02", 0, sr_b02);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b03", 0, sr_b03);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b04", 0, sr_b04);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b05", 0, sr_b05);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b06", 0, sr_b06);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_b07", 0, sr_b07);

                // Get quality Control
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_qc_500m", 0, sr_qc_500m);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_szen", 0, sr_szen);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_vzen", 0, sr_vzen);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_raz", 0, sr_raz);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_state_500m", 0, sr_state_500m);
                resHdf= getHDFMODPixelValueT (fMOD09A1Name, "sur_refl_day_of_year", 0, sr_day_of_year);

                // Extract combination data
                for (k=left_coo; k<=right_coo; k++)
                {
                    for (j=top_coo; j<=bot_coo; j++)
                    {
                        pos = k*dx + j;
                        data.sr_b01 = (int16)((long int)sr_b01[pos]);
                        data.sr_b02 = (int16)((long int)sr_b02[pos]);
                        data.sr_b03 = (int16)((long int)sr_b03[pos]);
                        data.sr_b04 = (int16)((long int)sr_b04[pos]);
                        data.sr_b05 = (int16)((long int)sr_b05[pos]);
                        data.sr_b06 = (int16)((long int)sr_b06[pos]);
                        data.sr_b07 = (int16)((long int)sr_b07[pos]);

                        data.sr_qc_500m = (uint32)((long int)sr_qc_500m[pos]);
                        data.sr_szen = (int16)((long int)sr_szen[pos]);
                        data.sr_vzen = (int16)((long int)sr_vzen[pos]);
                        data.sr_raz = (int16)((long int)sr_raz[pos]);
                        data.sr_state_500m = (uint16)((long int)sr_state_500m[pos]);
                        data.sr_day_of_year = (uint16)((long int)sr_day_of_year[pos]);

                        double nir = 1.0 * (data.sr_b02 * scale_factor_b02 + (*add_offset_b02) );
                        double swir = 1.0 * (data.sr_b06 * scale_factor_b06 + (*add_offset_b06) );
                        double blue = 1.0 * (data.sr_b03 * scale_factor_b03 + (*add_offset_b03) );
                        double red = 1.0 * (data.sr_b01 * scale_factor_b01 + (*add_offset_b01) );
                        
                        double ndvi = 1.0 * (nir - red) / (nir + red);
                        double lswi = 1.0 * (nir - swir) / (nir + swir);
                        double evi = 2.5 * (nir - red)/(nir + 6*red - 7.5*blue + 1);

                        /* Create query statement */
                        tmp_query.str("");
                        tmp_query << "INSERT INTO pixels (px_x, px_y, date_tkn ,sur_refl_b01, sur_refl_b02, sur_refl_b03, sur_refl_b04, sur_refl_b05, sur_refl_b06, sur_refl_b07, sur_refl_qc_500m, sur_refl_szen, sur_refl_vzen, sur_refl_raz, sur_refl_state_500m, sur_refl_day_of_year, ndvi, lswi, evi) VALUES (" << k << "," << j <<",'" << year << "." << date_taken <<"', " << data.sr_b01 << ", " << data.sr_b02 << ", " << data.sr_b03 << ", " << data.sr_b04 << ", " << data.sr_b05 << ", " << data.sr_b06 << ", " << data.sr_b07 << ", " << data.sr_qc_500m << ", " << data.sr_szen << ", " << data.sr_vzen << ", " << data.sr_raz << " , " << data.sr_state_500m << ", " << data.sr_day_of_year << ", " << ndvi << ", " << lswi << ", " << evi << ");";
                        query = tmp_query.str();
                        work W(C);
                        W.exec( query );
                        W.commit();

                    }
                    cout << "Image :\t" << year << "." << date_taken << "\t" << k << "/" << right_coo << " \n";
                }


                // free
                free (sr_b01);
                free (sr_b02);
                free (sr_b03);
                free (sr_b04);
                free (sr_b05);
                free (sr_b06);
                free (sr_b07);

                free (sr_qc_500m);
                free (sr_szen);
                free (sr_vzen);
                free (sr_raz);
                free (sr_state_500m);
                free (sr_day_of_year);

            }

        }
        C.disconnect ();
    }catch (const std::exception &e){
        cerr << e.what() << std::endl;
        return 1;      
    }

    ifs.close();

   return 0;
}





// ******************************************************
// Read HDF file

GByte getHDFMODPixelValueT(char *fMODISName, char *bName, int bNum, VOIDP outRaster){

    int32 SDidentifier; // Scientific Data identifier

    if (( SDidentifier = SDstart(fMODISName, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", fMODISName);
        exit(1);}

    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs);

    int out_index= -1;
    int i;

    int32 DataSetId; // Data Set identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    for(i=0; i < nDatasets; i++){
        if ((DataSetId = SDselect(SDidentifier, i)) == FAIL) {
            printf ("ERROR: Unables select dataset #%d\n", i);
            exit(2);}

        SDgetinfo(DataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes,
                  &DataSetDataType, &DataSetNumAttrs);

        if( !strcmp(DataSetName, bName)){
            out_index = DataSetId;
            //printf ("getHDFMODPixelValueT()--> DS: '%s'\n", DataSetName);
            i= nDatasets;}
        if (SDendaccess(DataSetId) == FAIL) {
            printf ("ERROR: Unable to close SDS\n");
            exit(3);}
    }   // end for

    if(out_index== -1){
        printf ("getHDFMODPixelValueT()--> ERROR: No Data Set named '%s'\n", bName);
        exit(4);}

    int32 edges[DataSetDimensionNumber]; // Array specifying the number of values to be read along each dimension
    int32 start[DataSetDimensionNumber]; // Array specifying from where start read values

    if (DataSetDimensionNumber == 3){
        start[0] = bNum;
        start[1] = 0;
        start[2] = 0;

        edges[0] = 1;
        edges[1] = DataSetDimensionSizes[1];    // 2030 in MOD021km
        edges[2] = DataSetDimensionSizes[2];    // 1354 in MOD021km
    }
    else {
        start[0] = 0;
        start[1] = 0;

        edges[0] = DataSetDimensionSizes[0];    // 2030 in MOD021km
        edges[1] = DataSetDimensionSizes[1];    // 1354 in MOD021km

    }

    intn SDcode;
    if((SDcode = SDreaddata(out_index, start, NULL, edges, (VOIDP)outRaster)) == FAIL){
        printf ("getHDFMODPixelValueT()--> ERROR: Loading error \n");
        exit(7);}

    return 0;
}

// ******************************************************
// Read Dataset Dimension

void findDimension(char *nameInputFile, int **arrayDimension, char *nameBand){

    int32 SDidentifier; // Scientific Data identifier
    if(( SDidentifier = SDstart(nameInputFile, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", nameInputFile);
        exit(1);}

    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs); // get information about the Input Dataset

    int32 DataSetId; // SubDataset identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    int i;
    int n;
    for(i=0; i< nDatasets; i++){ // access to all SubDatasets

        if((DataSetId = SDselect(SDidentifier, i)) == FAIL){ // if it's not possible open 'i' SubDataset
            printf ("ERROR: Unables select SubDataset #%d\n", i);
            exit(2);
        }
        SDgetinfo(DataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes, &DataSetDataType, &DataSetNumAttrs);
        if(0 == strcmp(DataSetName, nameBand)){
            *arrayDimension= (int *)malloc(sizeof(int)*DataSetDimensionNumber);
            for(n=0; n<DataSetDimensionNumber; n++)
                (*arrayDimension)[n]= DataSetDimensionSizes[n];
        }
    } // end Subdatasets for cycle

}

// ******************************************************
// Read Factor (scale and offset)


void getFactor(char * inputFileName, char * band, char *attribute, float **valueToSearch){

    int32 SDidentifier; // Scientific Data identifier

    if (( SDidentifier = SDstart(inputFileName, DFACC_READ)) == FAIL ){
        printf("ERROR: Unable to open inputfile '%s'\n", inputFileName);
        exit(1);}
    int32 nDatasets;
    int32 nGlobalAttrs;
    SDfileinfo(SDidentifier, &nDatasets, &nGlobalAttrs);

    int32 subDataSetId; // Data Set identifier
    char DataSetName[255]; // Name of the data set
    int32 DataSetDimensionNumber; // Number of dimensions in the data set; the maximum value is MAX_VAR_DIMS (or 32)
    int32 DataSetDimensionSizes[3]; // Array containing the size of each dimension in the data set; Y e X
    int32 DataSetDataType; // Data type for the data stored in the data set
    int32 DataSetNumAttrs; // Number of attributes for the data set

    int32 out_index= -1;
    int i;
    for(i=0; i < nDatasets; i++){
        if ((subDataSetId = SDselect(SDidentifier, i)) == FAIL) {
            printf ("ERROR: Unables select dataset #%d\n", i);
            exit(2);}

        SDgetinfo(subDataSetId, DataSetName, &DataSetDimensionNumber, DataSetDimensionSizes,
                  &DataSetDataType, &DataSetNumAttrs);

        if( !strcmp(DataSetName, band)){
            out_index = subDataSetId;
            // printf ("getReflectanceFactor()--> DS: '%s'\n", DataSetName);
            i= nDatasets;}
        if (SDendaccess(subDataSetId) == FAIL) { // to end a 'SDcreate' or a 'SDselect'
            printf ("ERROR: Unable to close SDS\n");
            exit(3);}
    } // end for

    if(out_index== -1){
        printf("getReflectanceFactor()--> ERROR: No Data Set named '%s'\n", band);
        exit(4);}

    int32 attr_index, num_type, count;
    char attr_name[64];
    int k;

    attr_index = SDfindattr(out_index, attribute);
    if (attr_index == FAIL){
        printf ("ERROR: SDfindattr failed %s\n", attribute);
        exit(5);}

    // printf("getReflectanceFactor()--> Found '%s' attribute\n", attribute);

    // Get information about the file attribute
    if (SDattrinfo(out_index, attr_index, attr_name, &num_type, &count) == FAIL){
        printf ("ERROR: SDattrinfo failed\n");
        exit(6);
    }

    *valueToSearch= (float *)malloc(sizeof(float) * count);
    if(valueToSearch== NULL){
        printf("ERROR: can not malloc()\n");
        exit(7);
    }

    float buffer[count];
    // --------------------------------

    // Read the attribute data
    if(SDreadattr(out_index, attr_index, buffer) == FAIL){
        printf ("ERROR: SDreadattr failed\n");
        exit(8);}
    
    for(k=0; k<count; k++){
        (*valueToSearch)[k]= (buffer[k]);
    }

    // Terminate access to the array
    int status = SDendaccess(out_index);
    if (status == FAIL){
        printf ("ERROR: SDendaccess failed\n");
        exit(9);}

    // Terminate access to the SD interface and close the file
    status = SDend(SDidentifier); // to end a 'SDstart'
    if (status == FAIL){
        printf ("ERROR: SDend failed\n");
        exit(10);}
}
// ******************************************************
